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Abstract 

In this article, we propose a Milstein finite difference scheme for 
a stochastic partial differential equation (SPDE) describing a large 
particle system. We show, by means of Fourier analysis, that the dis- 
cretisation on an unbounded domain is convergent of first order in 
the timestep and second order in the spatial grid size, and that the 
discretisation is stable with respect to boundary data. Numerical ex- 
periments clearly indicate that the same convergence order also holds 
for boundary-value problems. Multilevel path simulation, previously 
used for SDEs, is shown to give substantial complexity gains com- 
pared to a standard discretisation of the SPDE or direct simulation of 
the particle system. We derive complexity bounds and illustrate the 
results by an application to basket credit derivatives. 



1 Introduction 

Various stochastic partial differential equations (SPDEs) have emerged over 
the last two decades in different areas of mathematical finance. A classical 
example is the Heath- Jarrow- Morton interest rate model Heath et a/. (1992)] 
of the form 

d f . 1 ,t+x 2 



dr{x,t) = —\^rit,x) + -J^^'"a{t,u)dujdt + a{t,t + x)dWt,{^) 

where r{x,t) is the forward rate of tenor x at time t and a{t,t + x) its 
instantaneous volatility. In Benth fc Koekebakker(2008)| , a similar equation 
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has been proposed more recently to model electricity forwards. Most of 
the SPDEs studied share with ([T]) the property that the derivatives of the 
solution only appear in the drift term; in the case of ([T]) the volatility of the 
Brownian driver does not depend on the solution r at all. Numerical methods 
for hyperbolic SPDEs of the type ([T]) have been studied, for example, in 
|Roth(2002)| . 

This article, in contrast, considers the parabolic SPDE 

dv 1 d'^v dv 

where M is a standard Brownian motion, and fi and < p < 1 are real- valued 
parameters. It is clear that the behaviour of this equation is fundamentally 
different from those with additive or multiplicative noise. 

The significance of ([2]) for the following applications is that it describes 
the limiting density of a large system of exchangeable particles. Specifically, 
if we consider the system of SDEs 



dXl = fidt+y/l-p dWl + ^ dMt, (3) 

for 1 < i < A^, with (dWf, dW/) = 6ij and {dWi, dMt) = 0, where are 
assumed i.i.d. with finite second moment, the empirical measure 



i=l 



has a limit z/f for — > oo, whose density v satisfies ([2]) in a weak sense. For 
a derivation of this result in the more general context of quasi-linear PDEs 
Kurtz fc Xiong(1999)| . While the motivation in [Kurtz fc Xiong(1999)| 



see 



is to use a large particle system ([3]) to approximate the solution to the SPDE 
([2]), our view point is to use ([2]) as an approximate model for a large parti- 
cle system, and we will argue later the (computational) advantages of this 
approach in situations when the number of particles is large. 

As a first possible application, one may consider X* as the log price 
processes of a basket of equities, which have idiosyncratic components 
and share a common driver M (the "market"). If the size of the basket is 
large enough, the solution to the SPDE can be used to find the values of 
basket derivatives. In this paper, we study an application of a similar model 
to basket credit derivatives. 

We mention in passing that equations of the form (|2]) arise also in stochas- 
tic filtering. To be precise, ([2]) is the Zakai equation for the distribution of a 
signal X given observation of M, see e.g. [Bain fc Crisan(2009)| . 



2 



It is interesting to note that the solution to the SPDE ([2]) without bound- 
ary conditions can be written as the solution of the PDE 



du 1 , , d'^u du 

shifted by the current value of the Brownian driver, 

v{t,x) = u{t,x-^Mt). (5) 
In particular, if f(0,a;) = 5{xo — x), then 

viT X) - ' exp r (---o-^^T-^pMrr \ 

^'^"^Mi^'^H 2(i-p)T 

The intuitive interpretation of this result is that the independent Brow- 
nian motions have averaged into a deterministic diffusion in the infinite 
particle limit, whereas the common factor M, which moves all processes in 
parallel, shifts the whole profile (and also adds to the diffusion, via the Ito 
term) . 

In [Bush et a/. (2011)] , the analysis of the large particle system is extended 
to cases with absorption at the boundary (x = 0), 

Xi = 0, t>T^, 

T, = mf{t:Xi = 0}. (7) 

It is shown that there is still a limit measure Ut, which may now be decom- 
posed as 

ut = + Lt(5o, 

where 

1 ^ 

U = lim Lf = lim - ^ 1^.<, (8) 

i=l 

is the proportion of absorbed particles (the "loss function"), and the density 
V of satisfies ([2]) in (0, oo) with absorbing boundary condition 

v{t,0)=0. (9) 



Bush et g/. (201 1)1 consider applications to basket credit derivatives. For 



the market pricing examples, they consider a simplified model, where defaults 
are monitored only at a discrete set of dates. Between these times, the default 
barrier is inactive and ([2]) is solved on the real line by using (jl]) and ([5]). 
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For the initial-boundary value problem ([2]), however, such a semi- 
analytic solution strategy is no longer possible and an efficient numerical 
method is needed. Moreover, there is a loss of regularity at the bound- 
ary in this case, such that xuxx ^ L2 but not Uxx, as is documented in 
Krylov(1994)] . 



Recent papers on the numerical solution of SPDEs deal with cases rele- 
vant to ours, yet structurally crucially different. A comprehensive analysis 
of finite difference and finite element discretisations of the stochastic heat 
equation with multiplicative white noise and non-linear driving term is given 
in [Gyongy fc Nualart(1997)j|Gyongy(1999)| and |Walsh(2005)] , respectively. 



Lang(2010)| shows a Lax equivalence theorem for the SDE 

dXt = AXtdt + G{Xt)dMt, 

in a Hilbert space, driven by a process M from a class including Brownian 
motion, where A is a suitable (e.g. elliptic differential) operator and G a 
Lipschitz function. 

In this paper, we propose a Milstein finite difference discretisation for 
([2]) and analyse its stability and convergence in the mean-square sense by 
Fourier analysis. A main consideration of this paper is the computational 
complexity of the proposed methods, and we will demonstrate that a multi- 
level approach achieves a cost for the SPDE simulation no larger than that of 
direct Monte Carlo sampling from a known univariate distribution, 0(£:~^) 
for r.m.s. accuracy e, and is in that sense optimal. 

Multilevel Monte Carlo path simulation, first introduced in [Giles(2008a)| , 
is an efficient technique for computing expected values of path-dependent 
payoffs arising from the solution of SDKs. It is based on a multilevel de- 
composition of Brownian paths, similar to a Brownian Bridge construction. 
The complexity gain can be explained by the observation that the variance 
of high-level corrections - involving a large number of timesteps - is typically 
small, and consequently only a relatively small number Monte Carlo samples 
is required to estimate these contributions to an acceptable accuracy. Over- 
all, for SDKs, if a r.m.s. accuracy of e is required, the standard Monte Carlo 
method requires 0{e~^) operations, whereas the multilevel method based on 
the Milstein discretisation [Giles(2008b)] requires 0{e ^) operations. 

The ffist extension of the multilevel approach to SPDEs was for parabolic 
PDEs with a multiplicative noise term [Graubner ( 2008)] . There have also 



been recent extensions to elliptic PDEs with random coefficients Barth et a/. (2010) 
Cliffe et al.{20TT)\ . 



Our approach, for a rather different parabolic SPDE, is similar to the 
previous work on SDEs and SPDEs in that the solution is decomposed into 
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a hierarchy with increasing resolution in both time and space. Provided the 
variance of the multilevel corrections decreases at a sufficiently high rate 
as one moves to higher levels of refinement, the number of fine grid Monte 
Carlo simulations which is required is greatly reduced. Indeed, the total cost 
is only 0(e~^) to achieve a r.m.s. accuracy of e compared to an 0(£~^/^) cost 
for the standard approach which combines a finite difference discretisation 
of the spatial derivative terms and a Milstein discretisation of the stochastic 
integrals. 

The rest of the paper is structured as follows. Section [2] outlines the finite 
difference scheme used, and analyses its accuracy and stability in the stan- 
dard Monte Carlo approach. Section [3]presents the modification to multilevel 
path simulation of functionals of the solution. Numerical experiments for a 
CDO tranche pricing application are given in section HI providing empirical 
support for the postulated properties of the scheme and demonstrating the 
computational gains achieved. Section [5] discusses the benefits over standard 
Monte Carlo simulation of particle systems and outlines extensions. 

2 Discretisation and convergence analysis 

2.1 Milstein finite differences 

Integrating ([2]) over the time interval [t,t-|-A;] gives 

Making the approximation v{s, x) ~f (t, x) for t<s< t+k in the ffist integral 
and 

dv 

v{s,x)^v{t,x)~^— (Ms-Mt) 
in the second, and noting the standard Ito calculus result that 

{Ms-Mt) dM, = - ([^hTf-k) , 

where AM" = M^+fc - Mj = v^Z„ with ~ iV(0, 1), |Glasserman(2004) 



Kloeden fc Platen(1992)] , leads to the Milstein semi-discrete approximation 



. oil 1 / 

v-+\x)=v''{x)-{^ik+^Z^) — + -({l-p)k + pkZl 

ox / V 



dx"^ 
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Using a spatial grid with uniform spacing h, standard central difference 
approximations to the spatial derivatives Richtmyer &: Morton (1967)] then 
give the finite difference equation 



V 



+ ^^^^^|^e""«-2."+«;-o. (10) 

in which is an approximation to v{nk,jh). 

The spatial domain is truncated by introducing an upper boundary at 
Xmax = J h and using the boundary condition Wj = 0. Since the initial dis- 
tribution will be assumed localised, both the localisation error for a given 
path M, and the expected error of functionals of the solution, can be made 
as small as needed by a large enough choice x^ax > 0. 

The system of SDEs can then be written in matrix-vector form 



Hk + ^/^Zn {l-p)k + pkZl 
-2h 2h^ 



Vn+l = K - ^ .V -DlVn+ ' ^I^2K, (11) 



where Vn is the vector with elements f", j = 1, . . . , J — 1 and Di and D2 are 
the matrices corresponding to first and second central differences as explicitly 
given in Appendix lAl 

Remark 2.1 An alternative discretisation arises if the spatial derivatives 
are discretised first, and the Milstein scheme is applied to the resulting sys- 
tem of SDEs. The practical implication is that the ltd term then contains a 
second finite difference with twice the step size, D\ instead of D2, resulting in 
pentadiagonal discretisation matrices instead of tridiagonal ones, specifically 



Vn+i = Vn - - — ^ — ^^iK + ^^2v; + -^^^2 — -^i^n- (12) 

The accuracy and stability of the two schemes are similar, hence we will not 
go into details. 

To approximate the initial condition on the grid, we apply the initial 
measure to a basis of 'hat functions' (^j)o<j<j5 where 

\E'j(x) = — max {h — \x — Xj\, 0) , 

giving 

/■oo 

v'l = {"^j, Vo) = / ^j(x) Vo{x) dx. 
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For the corresponding PDE (p = 0), this is known to result in 0{h?) conver- 
gence provided k satisfies a certain stabihty hmit, even when the initial data 
are not smooth Pooley et a/. (2003), Carter fc Giles(2007)| . We will see that 
an extension of this analysis holds for SPDEs. 



2.2 Fourier stability analysis 

Finite difference Fourier stability analysis ignores the boundary conditions 
and considers Fourier modes of the form 

v^ = gneM^je), \e\<iT, (13) 

which satisfy ffTOj) provided 

gn+i = {ai9) + b{9) Z„ + c{9) Zl) g^, 



where 



T if^k 2{l-p)k . 2 

a(0] = 1 ; — smt^ — sm „, 

h 



biO) = ; — smb*, 

h 



Following the approach of mean-square stability analysis from Higham(2000) 
Saito fc Mitsui(1996)| , we obtain 



E[|^„+i|2] = E[{a + bZn + cZl){a* + b*Zn + c*Zl) \gn\^] 
= (|a+cp + |6p + 2|cp) E[kP], 

where a* denotes the complex conjugate of a. Mean-square stability therefore 
requires 



|ap + \b\^ + 3|cp + ac* + a*c 

1 - 4 sm^ I ^ ^ - (1 + 2/) ( A] ^ f " f (x) ^ + I ^ 



< 1, 



and enforcing this for all 6 leads to the two conditions summarised in the 
following theorem. 
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Theorem 2.1 The Milstein finite difference scheme ( fli]) is stable in the 
mean-square sense provided 

ff^k < 1-p, (14) 
^ < (l + 2p2)-^ (15) 

The analysis in Appendix |A] combines mean-square and matrix stabihty 
analysis to prove that in the limit k,h ^ the condition k/h"^ < (1 + 2p^)~^ 
is also a sufficient condition for mean-square stability of the initial-boundary 
value problem with the boundary conditions at j = and j = J. 

2.3 Fourier analysis of accuracy 

Fourier analysis can also be used to examine the accuracy of the finite dif- 
ference approximation in the absence of boundary conditions. Considering a 
Fourier mode of the form g{t) exp(mx), the PDE reveals that 

g{t) = g{0) exp(^-l{l-p)Kh-tK{pt+ ^pM^) , 

and therefore 

g{tn+l) = gifn) exp i(l-p)/t^/c - i K{pk+^/pk Zn)^ , 

where Mt^^-^ — Mt„ = ^JkZn- Fourier analysis of the discretisation gives 

5f„+i = {a{Hh) + h{H}i) Zn + c{Kh) Z^) g^, 
where aiff), b{9), c{6) are as defined before. Writing 

a{Kh) + h{Kh) Zn + c{Kh) Z"^ = exp ^— |(1 — p) t — i k {p k+^/fPk Zn) + e„j 

we obtain, after performing lengthy expansions using MATLAB's Symbolic 
Toolbox, 

Cn = -\/p^ K^pkZn - li \/f)k tC'Zn (1-p) k + pk Z"^- k^^ + rn, 
where, for all p > 1, 

E[|r„n < c{k,p) [e + khY ■ 
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When summing over T /k timesteps, ^ Z„ and ^ are both 0{k ' ) since 
they have zero expectation and 0{k~^) variance. Hence, it follows that 



E 



E 



1/2 



so the RMS error in the Fourier mode over the full simulation interval is 
0{k,h^). 

Following the method of analysis in Carter fc Giles(2007)| , it can be 
deduced from this that the RMS L2 and L^o errors for Vn are both 0{k, h?). 
This is consistent with the usual Oik^h?) accuracy of the forward-time central 
space discretisation of a parabolic PDE [Richtmyer fc Morton(1967)| , and 
the 0{k) strong accuracy of the Milstein discretisation of an SDE, see e.g. 
Kloeden fc Platen(1992)| . 



2.4 Convergence tests 

We now test the accuracy and stability of the scheme numerically. 

The chosen parameters for ([2]) are a = 0.22, p = 0.2, r = 0.042, /i = 
(r — cr^/2)/cr, v{0, x) =6{x—xo) with xq = 5. These values are typical for the 
applications later on. 

The upper boundary for the computation is x^ax = 16, and chosen to 
ensure that the use of zero Dirichlet data has negligible influence on the 
solution. 

We approximate the initial value problem without boundary conditions 
on [—16/3,16] (note = 5 is roughly in the centre), and also the initial- 
boundary value problem with zero Dirichlet conditions on [0, 16]. 

Figure [T] plots several solutions to the latter problem at T = 5, each 
corresponding to a different Brownian path Mt. It can be seen that for those 
realisations for which Mt has been largely positive, the boundary at x = 
has had negligible influence and so the solution is approximately equal to the 
displaced Normal distribution given by (jS]). 

For the unbounded case, we approximate the mean-square L2-error by 



E{h, kf 



E 



Y,{vf{uj)-v{Nk,3h-u)fh 



M J 



ll E - v{Nk,3h- Urn)? h, 



m=l j=0 



(16) 
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0.05 



Figure 1: Some solutions v{T,x) for different driving Brownian motions Wt- 



where J is the number of grid intervals, the number of timesteps, and the 
expectation is taken over Brownian paths u, of which Um are M samples. 

To study the convergence, we introduce decreasing grid sizes hi = ho 
and timesteps ki = ko4~'', motivated by the second order convergence in h 
and first order convergence in k as well as the stability limit for the explicit 
scheme, and denote Ei = E{hi, ki) the mean-square L2-error at level /. 

For the initial-boundary value problem, no analytical solution is known, 
but we can compute error indicators via the difference between a fine grid 
solution / with mesh parameters k and h, and a coarse solution c with mesh 
parameters 4k and 2h, 



e(h,ky 



E 



J/2 



M J/2 



(17) 



m=l j=0 



and define ei = e{hi, ki). 

On the coarsest level, ho = 4/3, /cq = 1/4, such that xq does not coincide 
with a grid point. 

Fig. [2] shows both the computed values of Ef and ef for the unbounded 
case, and for the bounded case. The results confirm the theoretical 
0{k,h'^) convergence in the unbounded case, and support the conjecture 
that the convergence order is unchanged in the bounded case. 
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Figure 2: Mean-square error estimators for the Milstein scheme. For the 
unbounded case, L2 difference between true and numerical solution, El as in 
(fT6|) . and coarse and fine grid solutions, ei as in ( IT7|) : for the bounded case, 
difference between coarse and fine grid solutions e/. 

We now formalise this conjecture about the error due to the SPDE dis- 
cretisation. Denote Ut = u{T, ■) the solution to the initial boundary value 
problem at time T for a given Brownian path and Ut its numerical approxi- 
mation with grid size h and timestep k (x h"^. 

Conjecture 2.1 The error in the solution at time T satisfies the strong error 
estimate in the L2 norm || • || 

y/E[\\UT-UT\n = 0{h'). 

Corollary 2.1 A Lipschitz payoff function P{Ut) has a similar strong error 
bound 

^n{p{UT)-P{UT)Y] = o{h'). 

We know Conjecture 12. II to be true for the initial value problem on M, and 
it conjectures that the introduction of the boundary condition at a: = does 
not affect the (weak and) strong error, which is supported by the numerical 
results. 

The conjecture assumes a;maa; = oo; in practice, a finite value for x^ax will 
introduce an additional truncation error which will decay exponentially in 
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Numerical tests with different parameters, not reproduced here, indicate 
that the error in the bounded case increases for values of Xq close to 0, 
in which case more paths contribute solutions with large higher derivatives 
close to the zero boundary, even if the asymptotic convergence order is still 
0{k, h"^). We analyse this loss of regularity further in Appendix [Bl 



3 Multilevel Monte Carlo simulation 

We now consider estimating the expectation of a scalar functional of the 
SPDE solution, an expected "payoff" P in a computational finance context. 

Defining P to be the payoff arising from a single SPDE approximation, the 
standard Monte Carlo approach is to average the payoff from independent 
SPDE simulations, each one using an independent vector of A^(0, 1) Normal 
variables Z = {Z^, Z^, . . . , Thus the estimator for the expected value 



IS 



M ^ 



N 

i=l 

The mean square error for this estimator can be expressed as the sum of two 
contributions, one due to the variance of the estimator and the other due to 
the error in its expectation. 



E 



Y - E[P] ■ ^ 



A^"^V[P] + (E[P]-E[P]'^ 



To achieve a r.m.s. error of e requires that both of these terms are 0(£:^). 
This in turn requires that = 0{e~'^), k = 0{e) and h = 0{e^^'^), based 
on the conjecture that the weak error E[P]— E[P] is 0{k,h'^). Since the 
computational cost is proportional to Nk~^h^^ this implies an overall cost 
which is 0(e-^/2)_ jj^g 

aim of the multilevel Monte Carlo simulation is to 
reduce this complexity to 0{e^'^). 

Consider Monte Carlo simulations with different levels of refinement, / = 
0, 1, . . . , L, with / = being the coarsest level, (i.e. the largest values for k 
and h) and level L being the finest level corresponding to that used by the 
standard Monte Carlo method. 

Let Pi denote an approximation to payoff P using a numerical discretisa- 
tion with parameters ki and hi. Because of the linearity of the expectation 
operator, it is clearly true that 

L 

E[Pi] =E[Po] + 5^E[Pi-Pi_i]. (18) 

1=1 
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This expresses the expectation on the finest level as being equal to the expec- 
tation on the coarsest level plus a sum of corrections which give the difference 
in expectation between simulations using different numbers of timesteps. The 
multilevel idea is to independently estimate each of the expectations on the 
right-hand side in a way which minimises the overall variance for a given 
computational cost. 

Let Yo be an estimator for E[Po] using A^o samples, and let for Z>0 be 
an estimator for K[Pi—Pi_i] using A^; samples. Each estimator is an average 
of Ni independent samples, which for l>0 is 

Ni 

Yi-N-'Y.iPl^'-P!^^- (19) 

i=l 

The key point here is that the quantity P/*'' — P/l\ comes from two discrete 
approximations using the same Brownian path. The variance of this sim- 
ple estimator is Y[Yi] = Nj~^Vi where Vi is the variance of a single sample. 
Combining this with independent estimators for each of the other levels, the 
variance of the combined estimator J2f=o ^ Xl^o -^z"^^- '^^^ correspond- 
ing computational cost is X]i=o where Q represents the cost of a single 
sample on level /. Treating the Ni as continuous variables, the variance is 
minimised for a fixed computational cost by choosing Ni to be proportional 
to \/Vi/Ci^ with the constant of proportionality chosen so that the overall 
variance is 0{e~'^). 

The total cost on level I is proportional to \/Vi Q. If the variance Vi 
decays more rapidly with level than the cost C/ increases, the dominant 
cost is on level 0. The number of samples on that level will be 0(e~^) and 
the cost savings compared to standard Monte Carlo will be approximately 
Cq/Cl, refiecting the different costs of samples on level compared to level 
L. On the other hand, if the variance Vi decays more slowly than the cost C; 
increases, the dominant cost will be on the finest level L, and the cost savings 
compared to standard Monte Carlo will be approximately Vl/Vq, reflecting 
the difference between the variance of the finest grid correction compared to 
the variance of the standard Monte Carlo estimator, which is similar to Vq. 

This outline analysis is made more precise in the following theorem: 

Theorem 3.1 Let P denote a functional of the solution of an SPDE for a 
given Brownian path M^, and let Pi denote the corresponding level I numerical 
approximation. 

If there exist independent estimators Yi based on Ni Monte Carlo samples, 
and positive constants a, ^, 7, ci, C2, C3 such that a > | 7 and 
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i) E[Pz-P] 



< ci2 



-al 



E[Po] 



/ = 



III) N[Yi] < C2N-^2-P' 

iv) Ci < c^Ni2'^\ where Ci is the computational complexity ofYi 

then there exists a positive constant C4 such that for any e < e^^ there are 
values L and Ni for which the multilevel estimator 



/=0 



has a mean-square-error with bound 



MSE = E 



Y - E[P] 



< e 



with a computational complexity C with bound 

CaE-'^, /3>7, 
C4£-2(loge)2, /3 = 7, 
^^^-2-(r-/3)/«^ 0</3<7. 



C < 



Proof The proof is a slight generahsation of the proof in [Giles(2008a)] . □ 

In our application, we choose hi oc and ki oc so that the ratio 
ki/hf is held fixed to satisfy the finite difference stability condition. The 
computational cost increases by factor 8 in moving from level I to / + !, so 
7 = 3. 

Given that the payoff is a Lipschitz function of the loss approximations 
at various dates, Conjecture 12. II implies that the weak error is also OQi^) and 
so a = 2 > 1 7. Also, due to the triangle inequality 



V[P,-P;_i] < ^V[Pi-P] + ^V[Pz_i-P] 

< v/E[(Pi-P)2] + A/E[(P;_i-P)^ 



Conjecture 12. II and its corollary give 13 = A. Consequently, the computational 
cost to achieve a r.m.s. error of e is 0(e~^). 
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4 Numerical experiments 



In this section we study the numerical performance of the algorithms pre- 
sented earlier on the example of CDO tranche pricing in a large basket limit. 

4.1 CDO pricing in the structural credit model 

Basket credit derivatives provide protection against the default of a certain 
segment ('tranche') of a basket of firms. A typical example is that of a 
collateralized debt obligation where the protection buyer receives a notional 
amount, minus some recovery proportion < i? < 1, if firms in a specified 
tranche of the basket default, and in return pays a regular spread until the 
default event occurs. 

The arbitrage-free spread depends crucially on the (risk-neutral, when 
hedged with defaultable bonds) probability of joint defaults. For a tranche 
with attachment point < a < 1 and detachment point 1 > d > a, the 
outstanding tranche notional 

P{Lt) = max(ci - Lt, 0) - max(a - Lt, 0), (20) 

where the loss variable Lt is the proportion of losses in the basket up to time 
t, determines the spread and default payments related to that tranche. The 
risk-neutral value of the tranche spread can be derived as 

e-^-E^[P(L^^_J-P(L^J] 



see, e.g., |Bush et al.{20TT)\ . Here, n is the maximum number of spread 



payments, T the expiry, Tj the payment dates for I < i < n, 5 = 0.25 the 
interval between payments. Spreads are quoted as an annual payment, as a 
ratio of the notional, but assumed to be paid quarterly. 

For an extensive survey of credit derivatives and pricing models we refer 
the reader to |Schonbucher(2003)| , and note only that they typically fall in 
one of two classes: so-called reduced-form models, which model default times 
of firms directly as (dependent) random variables; structural models, which 
capture the evolution of the firms' values, and model default events as the 
first passage of a lower default barrier. We will consider the latter here. 

In a structural model in the spirit of the classical works of |Merton(1974)| 
and Black fc Cox(1976)] , as extended to multiple firms e.g. in the work of 



^ There is a variation for the equity tranche, and recently sometimes the mezzanine 
tranche, but we do not go into details here. 



15 



Hull et a/.(2005)] , a company z's firm value, 1 < i < A^, is assumed to follow 



a model of the type 
d^ 



A] 



dt + ^/l-pi (Ti dW^ ^^iO, dMt, = a 







where < pi < 1 is a correlation parameter, r is the risk-free interest rate, 
Oi the volatility, and M, are standard Brownian motions. 

Here, the individual firms are correlated through a common 'market' fac- 
tor M, but independent conditionally on M. We make in the following the 
assumption that the firms are exchangeable in the sense that their dynamics 
is governed by the same set of parameters, specifically p = Pi, a = ai. Their 
initial values are not necessarily identical which allows for different default 
probabilities for individual firms, consistent with their CDS spreads. 

In this framework, the default time Tq of the i-th firm is modelled as 
the first hitting time of a default barrier 5*, for simplicity constant, and the 
distance-to-default 

X; = -(logA^-logS^), (22) 
a 

evolves according to ([3]), where p = (r — |cr^) /a, x* = (log a* — logS*) /a. 
Tq as in ([7]) is precisely the default time. 

In the majority of applications, Monte Carlo simulation of the firm value 
processes is used for the estimation of tranche spreads, see e.g. Hull et a/.(2005)| 
or 



Fouque et a/. (2008), Carmona et a/.(2009)| . This is largely enforced by 



the size of iV, for instance in the case of index tranches = 125. 

This, however, puts the model precisely in the realm of the large basket 
approximation ([2]). See Bujok k. Reisinger(2012)| for a numerical study of 



this large basket approximation. The loss functional is thereby approximated 
by 

in terms of the numerical SPDE solution v^'^^ at time Tj, which feeds into 
the estimator for the outstanding tranche notional fl20l) and subsequently the 
tranche spreads 



4.2 Pricing results 

All following results are for a representative set of parameter values, taken 
from a calibration performed in [Bush et a/.(201lj] to 2007 market data, 
a = 0.22, p = 0.2, r = 0.042, p = (r-a^/2)/a. While the initial distribution 
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used in Bush et a/. (2011)] is somewhat spread out to match individual CDS 



spreads of obhgors, most of the mass is around Xq = 5 and for simphcity we 
centre all firms there for the numerical tests, v{0,x) = 6{x — xo). The upper 
boundary is chosen as Xmax = 16. 

We consider a maturity T = 5 and tranches [a,d] = [0,0.03], [0.03,0.06], 
[0.06,0.09], [0.09,0.12], [0.12,0.22], [0.22,1]. 

Figure [3] shows the multilevel results for the expected protection payment 
from the first tranche, 

n 

5^e-'^^»E^[P(LT._J-P(LTj], 

i=l 

expressed as a fraction of the initial tranche notional. 

We pick mesh sizes hi = /io2~' and ki = fco4~' for Z > and ho = 8/5, 
ko = 1/4. This is motivated by the second and first order consistency in h 
and k respectively, and within the stability limit f lTSj) . 

The top-left plot shows the convergence of the variance Y[Pi — Pi-i] as 
well as the variance of the standard single level estimate, while the top-right 
plot shows the convergence of the expectation E[P/ — P/_i]. 

The bottom-left plot shows results for different multilevel calculations, 
each with a different user-specified accuracy requirement. For each value of 
the accuracy e, the multilevel algorithm determines (by formula (12) from 
|Giles(2008a)| ) the numbers of levels of refinement which are needed to ensure 
that the contribution to the Mean Square Error (MSE) due to the weak 
error on the finest grid is less than and it determines (by formula (10) 
from |Giles(2008a)| ) the optimal number of samples on each level so that the 
combined variance of the multilevel estimate is also less than ^e^, and hence 
the MSE is less than e^. 

Treating a single finite difference calculation on the coarsest level as hav- 
ing a unit cost, the bottom-right plot compares the total cost of the standard 
and multilevel algorithms. Since the objective is to achieve a computational 
cost which is approximately proportional to e~'^, it is the cost C multiplied by 
which is plotted versus e. The results confirm that e'^C does not vary much 
as £ ^ for the multilevel method, in line with the prediction that e'^C oc 1, 
whereas it grows significantly for the standard method since e'^C oc 8^^, where 
L is the index of the finest level. Hence, there is a big jump in the cost of 
the standard method each time it is necessary to switch to a finer level to 
ensure the weak error is less than e/v^, whereas the jump is minimal for the 
multilevel method. 

In practice, default is not monitored continuously but only at a discrete 
set of times, and for pricing often the simplifying assumption is made that 
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Figure 3: Multilevel results for the expected loss from the first tranche. 
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default is only determined at the spread payment dates Ti. 

This can be incorporated as follows. In the time intervals (T/_i,Ti), we 
solve ([2]) on a sufficiently large domain (e.g. [—4, 16]), and apply the following 
interface conditions at Tf. 

Ymvvit x) = \ 0, X < 0, -^gN 

t\.Ti ^' \\mvt^TiV{t,x), X > 0. 

To maintain quadratic grid convergence in spite of the discontinuity in- 
troduced by (125]) . we choose the mesh such that a grid point coincides with 
the boundary (e.g. by setting /iq = 2 in the above example), and set the 
numerical solution after default monitoring to for grid coordinates below 
and to 1/2 its previous value at 0, see e.g. [Pooley et a/.(2003)| . 



It is seen in Fig. H] that convergence is very similar to the previous case. 



5 Conclusions 

5.1 Complexity and cost 

Here we discuss the computational complexity of the multilevel solution of 
the SPDE compared with the alternative use of the multilevel method for 
solving the SDEs which arise from directly simulating a large number of 
SDEs. 

We have already explained that to achieve a r.m.s. accuracy of e requires 
0{e~'^/'^) work when solving the SPDE by the standard Monte Carlo method, 
but the cost is 0{e~'^) using the multilevel method, provided Conjecture 12.11 
is correct. 

Consider now the alternative of using a finite number of particles (firms), 
M, to estimate the tranche loss in the limit of an infinite number of particles 
(firms). In this case, empirical results suggest that there is an additional 
0{M~^) error (see also Bujok fc Reisinger(2012)| ), and the proof of this 



convergence order is the subject of current research. Taking this to be the 
case, the optimal choice of M to minimise the computational complexity 
to achieve an r.m.s. error of e is 0{e~^). Using the standard Monte Carlo 
method, the optimal timestep is 0{e), and the optimal number of paths is 
0(5"^), so the overall cost is 0{e~'^). Using the multilevel method for the 
SDEs reduces the cost per company to 0(6"^), so the total cost is 0{e~^). 

This complexity information is summarised in Table [TJ There is also a 
practical implementation aspect to note. The computational cost per grid 
point in the finite difference approximation of the SPDE is minimal, requiring 
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accuracy e 



Figure 4: Multilevel results for the expected loss from the first tranche, with 
discrete monitoring. 
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Table 1: Comparison of the complexity of the SDE and SPDE models using 
the standard and multilevel Monte Carlo methods 



Method / model 


SDE 


SPDE 


Standard MC 






Multilevel MC 







just three floating point multiply-add operations if equation flTU]) is re-cast 
as 

v^^' = at;;_i + bvj + cvj^, 

with the coefficients a, b, c computed once per timestep, for all j. 

If we let C be the cost of generating all of the Gaussian random num- 
bers Zn for a single SPDE simulation, then the cost of the rest of the finite 
difference calculation with 20 points in x (as used on the coarsest level of 
our multilevel calculations) is probably similar, giving a total cost of 2C for 
each SPDE. On the other hand, each SDE needs its own Gaussian random 
numbers for the idiosyncratic risk, and so the cost of simulating each SDE is 
approximately C, roughly half of the cost of the SPDEs on the coarsest level 
of approximation, giving a total cost of M C . 



5.2 Further work 

We have shown that stochastic finite differences combined with a multilevel 
simulation approach achieve optimal complexity for the computation of ex- 
pected payoffs of an SPDE model. In the case of an absorbing boundary, the 
complexity estimate is a conjecture in so far it relies on the convergence order 
of the finite difference scheme, which does not follow from the Fourier analy- 
sis of the unbounded case. The matrix stability analysis in Appendix \K\ could 
form part of a rigorous analysis if a Lax equivalence theorem could be proved. 
In the case of multiplicative white noise this is shown in Lang (2010)| . One 



difficulty in the present case of an SPDE with stochastic drift is the loss of 
regularity towards the boundary, which may be accounted for by weighted 
Sobolev norms of the solution, but even then it is not clear that convergence 
of the functionals of interest follows. 

There are several possible extensions of the present basic model as dis- 
cussed in Bujok fc Reisinger(2012)| , ranging from stochastic volatihty and 



jump-diffusion to contagion models. The methods developed in this paper 
should be of use there also, building for example on multilevel versions of 
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jump-adapted discretisations for jump-diffusion SDEs Bruti-Liberati & Platen(2010) 
Xia fc Giles(20T2)] . 
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A Mean-square matrix stability analysis 



If Vn is the vector with elements Vj',j = 1, . . . , J— 1 then the finite difference 
equation can be expressed as 



K+i 
A 

B 
C 



{A+BZn + CZD Vn, 

uk ^ (l—p) k ^ 



/f)k 
~2h' 



pk 



where / is the identity matrix and Di and D2 are the matrices corresponding 
to central first and second differences, which for J = 6 are 



/ 1 \ 

-1 1 

-10 1 

-1 1 

V -10/ 



/ -2 1 \ 

1 -2 1 

1 -2 1 

1 -2 1 

V 1-2/ 



Prom the recurrence relation we get 

= ElV^iA-^ + B^Zn + C^ZDiA + BZn + CZDVn] 
= E [Vf {{A+Cf{A+C) + B'^B + 2 C^C ) K] ■ 

Noting that Di is anti-symmetric and D2 is symmetric, and that 

D1D2 - D2D1 = Ei- E2, Df = D3 + Ei + E2, 

where D3 corresponds to a central second difference with twice the usual 
span, 

/ -3 1 \ 
0-201 
^3= 10-201 
10-20 
V 10-3/ 

(with the end vahics of —3 being chosen to correspond to V^i = —Vi and 
Vj+i = — Vj_i), and Ei and E2 are each entirely zero apart from one corner 
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element, 



/2 



El 



\ 



Eo 



\ J 

then after some lengthy algebra we get 



v 



2/ 



E [If ((A+C)^(A+C) + B^B + 2 C^C ) K] 
= E [Vf MK] - (ei + 62) E[K)2] - (ei - e^) E[K_i; 



where 



and 



M = I - —D2 + 



Ah* 



pk p^k'^\ 

+ -rrrr I ^3, 



ei 



.2u2 
2h^ ' 2/^2"' 



4/i2 4/;.2 



62 



2/i3- 



It can be verified that the m^^ eigenvector of M has elements sinj^^ for 
9m = mn/J, and the associated eigenvalue is 

\a{em) + c{em)\^ + + 2|c(^™)|2, 

where a{6) , b{6) , c{6) are the same functions as defined in the mean-square 
Fourier analysis. In addition, in the limit h,k/h — )■ 0, ei±e2 > 0, and 
therefore in this limit the Fourier stability condition 



sup [\a{9) + c 



2|c 



is also a sufficient condition for mean-square matrix stability. 



B Regularity considerations 

Figure [5] shows the convergence behaviour as the computational grid is re- 
fined. Level has /i= 1/4, = T/4; h is reduced by factor 2 and k by factor 
4 in moving to finer levels. 

The top-left plot shows the convergence of E[P/ — P^-i], with notation as 
in Sections [3] and the top-right plot shows the convergence of V[P/ — P/_i]. 

The bottom two plots show the behaviour of d^v/dx^{T,0), which is 
estimated on each grid level using the standard second difference 

d'^v , , _ V2 - 2f 1 + vq 
dx^ ^ ' ^ h? 
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Figure 5: Convergence plots as a function of grid level 

The left plot indicates that the mean of this quantity is well behaved, but 
the right plot indicates a singular behaviour of its variance, with the value 
increasing rapidly with increased grid resolution. 

This is in accordance with the result shown in Krylov(1994)| , that for 



the unique solution v to the SPDE, x v^x is square integrable, but v^x has a 
singularity at 0. 
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